      REAL*8 FUNCTION GHANGL(TSEC50,KPRECES)
      IMPLICIT REAL*8(A-H,O-Z)
C
C  THIS ROUTINE COMPUTES THE GREENWICH HOUR ANGLE AT A GIVEN TIME FOR
C  EITHER A FIXED OR PRECESSING EQUINOX.
C
C  VARIABLE DIM  TYPE  I/O DESCRIPTION
C  -------- ---  ----  --- -----------
C
C  TSEC50    1   R*8    I  THE TIME FOR WHICH THE GREENWICH HOUR ANGLE 
C                          IS WANTED. UNITS ARE MEAN SOLAR SECONDS 
C                          SINCE 1/1/50, 0.0 HRS UT.
C                          TSEC50 MAY BE ANY VALUE, + OR -.
C
C                          TSEC50 SHOULD BE COMPUTED USING THE JULIAN
C                          DATE(UT) OF INTEREST MINUS THE JULIAN DATE
C                          OF 1/1/50 0.0 (=2433282.5) MULTIPLIED BY
C                          86400. THAT IS, DO NOT ADD IN LEAP SECONDS.
C
C  KPRECES   1   I*4    I  A FLAG INDICATING THE COORDINATE SYSTEM 
C                          TO WHICH THE HOUR ANGLE IS REFERENCED.
C                          = 0, MEAN EQUATOR AND EQUINOX OF 1950.0
C                          = OTHERWISE, MEAN OF DATE.
C
C  GHANGL    1   R*8    O  THE GREENWICH HOUR ANGLE MEASURED IN A 
C                          POSITIVE SENSE(EASTWARD) IN THE EQUATORIAL 
C                          PLANE FROM THE EQUINOX(X-AXIS) TO THE 
C                          GREENWICH MERIDIAN. IN RADIANS.
C
C**********************************************************************
C
C  CODED BY C PETRUZZO, 11/82.
C     MODIFIED....  7/83. CJP. COMMENT CLARIFICATION. CODE CLEANUP, AND
C                         CORRECTION TO METHOD OF HANDLING THE
C                         NON-PRECESSING EQUINOX.
C
C**********************************************************************
C
C
C  METHOD: AS IN THE ASTRONOMICAL ALMANAC, 1983, PAGE B-6. THOSE 
C          FORMULAE ARE FOR MEAN OF DATE.
C
      REAL*8 SECDAY/86400.D0/,SECHR/3600.D0/,MINHR/60.D0/
      REAL*8 TWOPI/ 6.283185307179586D0 /
      REAL*8 DEGRAD/ 57.29577951308232D0 /
      INTEGER INIT/1/
C
C
      IF(INIT.NE.0) THEN    ! INITIALIZATION ON FIRST CALL.
        INIT=0
        CENT50=0.5D0  ! CENTURIES 1/0/1900, 12 HR, TO 1/1/50, 0 HR (UT).
        RADHR = 15.D0/DEGRAD   ! 15 (DEG/HR) / (DEG/RAD)= RAD/HR.
C       SET COEFFICIENTS FOR RU FORMULA, MEAN OF DATE CASE.
        C1=18.D0 + 38.D0/MINHR + 45.836D0/SECHR
        C2=8640184.542D0/SECHR
        C3=0.0929D0/SECHR
        END IF
C
C    CENT = FRAC OF A JULIAN CENTURY(MEAN SOLAR DAYS) FROM 1/0/00, 12 HR
      CENT = (TSEC50/SECDAY)/36525.D0 + CENT50  
C    UT = UNIVERSAL TIME. HOURS.
      UT = DMOD(TSEC50/SECHR,24.D0)  
C    RU = RACN OF FICTITIOUS MEAN SUN. HOURS.
      RU = C1 + C2*CENT + C3*CENT*CENT
C    GMST = GREENWICH MEAN SIDEREAL TIME(GHA OF MEAN EQNX OF DATE). HRS.
      GMST = UT + 12.D0 + RU
C    GHA = GREENWICH HOUR ANGLE. RADIANS.
      GHA = GMST * RADHR
C
C    GHA IS NOW IN MEAN OF DATE. CONVERT TO MEAN OF 1950.0 IF NEEDED.
      IF(KPRECES.EQ.0) THEN
        SHIFT=PRECES50(TSEC50,1)
      ELSE
        SHIFT=0.D0
        END IF
      GHA = GHA - SHIFT
C
C  WRAP IT UP.
      GHANGL = DMOD(GHA,TWOPI)
C
      RETURN
      END
